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Abstract 

We present an alternative derivation of the dynamical density functional theory for the one 
body density profile of a classical fluid developed by Marconi and Tarazona [J. Chem. Phys., 110, 
8032 (1999)]. Our derivation elucidates further some of the physical assumptions inherent in the 
theory and shows that it is not restricted to fluids composed of particles interacting solely via 
pair potentials; rather it applies to general, multi-body interactions. The starting point for our 
derivation is the Smoluchowski equation and the theory is therefore one for Brownian particles 
and as such is applicable to colloidal fluids. In the second part of this paper we use the dynamical 
density functional theory to derive a theory for spinodal decomposition that is applicable at both 
early and intermediate times. For early stages of spinodal decomposition our non-linear theory 
is equivalent to the (generalised) linear Cahn-Hilliard theory, but for later times it incorporates 
coupling between different Fourier components of the density fluctuations (modes) and therefore 
goes beyond Cahn-Hilliard theory. We describe the results of calculations for a model (Yukawa) 
fluid which show that the coupling leads to the growth of a second maximum in the density 
fluctuations, at a wavenumber larger than that of the main peak. 

PACS numbers: 05.20.Jj, 64.60.My, 61.20.Lc 
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I. INTRODUCTION 



Classical density functional theory (DFT)A is a remarkably successful theory for describing 
the rich behaviour of the equilibrium structure and thermodynamics of fluids in external 
potentials. DFT has been used to describe a wide variety of fluid interfacial, confinement 
and even freezing phenomena; for example, DFT has been a vital tool in understanding the 
wetting behaviour and surface phase transitions of fluids adsorbed on various substrates^. 
Given the success of DFT for describing static inhomogeneous fluids, it is of great interest 
to be able to build upon and incorporate these theories into a theory for the dynamics of 
inhomogeneous fluids. 

There have been several approaches to obtaining an equation of motion for the one 
body density profile p(r, t) of a classical fluid. The form that these theories takes depends 
somewhat on how the particular authors define p(r,t). For some, p(r,t) is an "ensemble" 
average over the possible configurations of the system at time t, given an ensemble of starting 
configurations at an earlier time t = 0. Using this definition for p(r,t), there is clearly a 
unique density profile p(r, t) at a given time t, and therefore the equation governing the 
dynamics of this density profile will be deterministic. This is the philosophy behind the 
approach of Marconi and Tarazona (MT) in Refs.^. Their approach obtains more formally 
some of the results proposed by earlier authors, such as Evans^ and Dieterich et a/.—, where 
it is assumed from the outset that the gradient of the chemical potential V/x(r, t) is the 
thermodynamic force driving a particle current 

j(r,t) = -r P (r,t)VMr,t), (1) 

where T is a mobility constant. An expression for fi is obtained within DFT by assuming 
that, as in the case of the equilibrium fluid, the chemical potential is given by the functional 
derivative of the Helmholtz free energy functional with respect to the density profile^. This 
assumption, together with Eq. and the continuity equation 

dp(r,t) 



Of 



-V.j(r,t), (2) 



provide a basis for the deterministic formulation of dynamical DFT (DDFT). Equations of 
this form have been used, for example, to study the dynamics of freezing^ and of solvation^. 
More recently, the more systematic approach of MT has been used with much success to 
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describe the dynamics for several different systems. These applications refer to particles in 
various external, time dependent, potential o 9 ' 10 ' 11 ' 12 . For the systems considered, the agree- 
ment between theory and the results from Brownian dynamics simulations have generally 
been very good. 

An alternative approach to obtaining a DDFT is to view the fluid one body density pro- 
file (denoted p(r,t) in order to distinguish it from p(r, t), the "ensemble" averaged density 
profile) as some sort of spatial and/or time coarse grained average of the density operator 
p(r,t) = ~52i=i <K r — r i(t))i where the rj(i) are the positions of the iV particles in the sys- 
tem. In this case, for Brownian particles, the equation governing the dynamics of p(r, t) 
will, of course, still contain a stochastic element. This is the viewpoint of a number of 
theorie o 13 ! 14 ! 15 ! 16 ! 17 ! 18 . There is some confusion in the literatur e - 18 ! 19 (see also Ref.—) as to 
what precisely is meant by p(r,t). We will attempt to clarify these issues elsewhere 21 . 

In this paper we adopt the deterministic viewpoint of MT and consider p(r, t) to be an 
ensemble average, i.e. an average over all realisations of the stochastic noise in the preceding 
interval, until time t. We employ the Smoluchowski equation as the starting point for an 
alternative derivation of the DDFT for classical fluids developed by MT in Refs.^. The 
present scheme for deriving the DDFT bears some similarity in spirit to that given recently 
in Ref.—, where projector-operator techniques are used to obtain first the Smoluchowski 
equation and subsequently an equation of motion for the fluid one body density profile - 
the DDFT. Before proceeding with our derivation of the DDFT, in Sec. |H]we give a brief 
introduction to the Smoluchowski equation, expounding some of the physical assumptions 
concerning the dynamics of the fluid that are implicit in this equation. In Sec. IIHI we 
proceed with the derivation of the DDFT of MT from the Smoluchowski equation. Sec. IIVI 
describes an important application of the DDFT to analyse the short and intermediate-time 
dynamics of spinodal decomposition relevant to colloidal fluids. Finally in Sec. we make 
some concluding remarks. 

II. THE SMOLUCHOWSKI EQUATION 

The Smoluchowski equatio n 23 ! 24 ' 25 ! 26 ' 27 is a Fokker-Planck equation (or generalised diffu- 
sion equation) for interacting Brownian particles. A physically intuitive way of arriving at 
this equation, e.g. Ref.—, proceeds as follows: For a fluid of N Brownian particles, one imag- 
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ines applying a force on the particles, where the force on the j th particle is Fj = —Vj^(r N , t), 
(r N = {ri, r2...rAr} is the set of position coordinates for the N particles), so that the system 
is prevented from relaxing to its equilibrium distribution. The non-equilibrium probability 
density function in this situation will be: 

P(r N , t) = i exph/^r*, t) - /3U(r N , t)} (3) 

where Z is a normalisation factor, (3 = 1/ksT is the inverse temperature and U(r N ,t) is the 
potential energy due to the interparticle interactions and any external potentials. Taking 
the gradient of Eq. (jSJ) we obtain: 

F s = VMT»,t)+k B T ?0^- . (4) 

If Fj is switched off, there will be a force — Fj driving the diffusion of particle j. We now 
assume that for time scales ^> rg, the Brownian time scale, the velocity of the i th particle 
is 

N 

(5) 

3=1 

where = /3D^(r , t) is the mobility tensor and is the diffusion tensor. It is implicitly 
assumed that as the particles interact, the momentum degrees of freedom equilibrate much 
faster than the positional degrees of freedom, and we have effectively averaged over the 
momentum degrees of freedom, whilst keeping the particle coordinates fixed. For a colloidal 
fluid, this thermal equilibration should occur via the solvent, and this approximation should 
be a good one to make. For atomic fluids, this may not be the case, especially for particles 
interacting via harshly repulsive (hard-sphere like) potentials. For fluids interacting with 
softer potentials, such as in the Gaussian core model^, this might be a reasonable approxi- 
mation to make, particularly when the fluid is at high densities where each particle interacts 
with a large number of neighbours and so the momentum degrees of freedom can equilibrate 
faster. Since the particle number is conserved, we can expect the fluid to obey the continuity 
equation: 

^^ = -f:V,[v l P(r-,t)]. (6) 

i=l 
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Substituting Eqs. (HJ) and (J3J) into one finds 



dP{T»,t) A A 

i=l j=l 



dt 



££v,iv 



k B TV 



j 



+ V J U{v\t) P(r N ,t). (7) 



If the potential energy term U(r N ,t) = 0, then = FS^ = f3D5ij, where D is 
the diffusion coefficient and Eq. (J7j) reduces to the diffusion equation, (d/dt)P(r N , t) = 
f3T VfP(r Ar , t). For a system of interacting particles the diffusion tensor does not, in 
general, take such a simple form 2 ^ However, if we neglect the hydrodynamic interactions, 
we can replace T^- by its 'mean- field' value, T5ij, and then Eq. (JJJ) reduces to a generalised 
diffusion equation, termed the Smoluchowski equation: 

d P( !L V) = rf>,.[^TV, + V^(rV)]P(rV). (8) 



i=l 



More formally, the Smoluchowski equation is the Fokker-Planck equation for a system of N 
Brownian particles in the large friction limi t 23 i 26 i 27 . The Langevin equation for a system of 
N Brownian particles of mass m is: 

+ T- 1 ^ = -V^(rV) + Mt), (9) 

where Xj(i) = £f (£), £f(£)) is a white noise term with the property (£"(£)) = and 

= 2k B T5ij5 au 5(t — t'). When the friction constant T -1 is large, we may neglect 
the second derivative with respect to time in Eq. 0, and we obtain the stochastic equation 
of motion^: 

^ = -rV<l/(r* t) +VX i (t). (10) 
The (generalised) Fokker-Planck equation for the distribution function P(r , t) correspond- 



ing to this Langevin equation is Eq. 



k 23.26.27 



III. DYNAMICS OF THE ONE-BODY DENSITY PROFILE 

In this section we derive an equation for the time evolution of the one-body density profile, 
p(r,t), from the Smoluchowski equation, Eq. (|8*|) . For a similar approach based solely on 
pair potentials see Refs. 29 ' 30 . The one body density is merely the integral of the probability 
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distribution function: 



p(ri,t) = N J dr 2 ... J dr N P(r N ,t) (11) 
Similarly, the two-body density is 

p( 2 >(n,r 2 ,t) = N(N-l) J dr 3 ...y dr N P(r N ,t), (12) 

and in general the n-particle density is 

P (n VV) = ^^yy J dv n+1 ...J dr N P(r N ,t). (13) 



Using Eqs. (jll)) - (jT3j) . and assuming that the potential energy function can be expressed 
in terms of a one-body external potential acting on each particle, V ex t(ri,t), and that the 
particle interactions are a sum of pair potentials, t> 2 (rj, i\,-), three body potentials f 3 (i"j, rj, r^), 
and higher body interactions: 

N 1 N 

u(v N ,t) = J2 v ^i,t) + ^EE^^) 

i=l j^i i=l 

1 N 

+ gE 2222v 3 {T it T j} T h ) + ... (14) 

then we find that on integrating Eq. (JBJ), one obtains 

r 7^ = kuTV-LpiTut) 

+ Vi.|>(r 1 ,t)ViV r B Bt(ri,t)] 
+ Vi. / dr 2 p (2) (n, r 2 , t)Vi^ 2 (r 1 , r 2 ) 

+ Vi. J dr 2 y dr 3 p (3) (ri,r 2 ,r 3 ,t)Vit;3(r 1 ,r 2 ,r 3 ) 

+ .... (15) 

We note that if dp(r,t)/dt = 0, then Eq. (|T5j) is equivalent to the derivative of the first 
equation of the YBG hierarchy^. 

For a fluid in equilibrium, there is an exact sum rule^ which relates the gradient of the 
one body direct correlation function to the interparticle forces acting on a particle (recall 
that — fc^Tc^^r) is the effective one-body potential due to interactions in the fluid). If the 
particles interact solely via pair potentials: 

-k B Tp(r)Vc {1 \r) = J dr'p (2 )(r, r')Vti 2 (r, r'). (16) 



This result can be generalised straightforwardly to fluids where the particles interact via 
many-body potentials, Eq. (JHJ), giving: 



fc B Tp(n)VcW(ri) = 

J2 / dr 2 ... / dr n p< n >(r n )Vi« n (r n ). (17) 

n=2 J J 



From equilibrium statistical mechanics one also knows that (r) is equal to the functional 
derivative of the excess (over ideal) part of the Helmholtz free energy functional 5 : 

c(1)(r) = _ ;3 ^)] i (18) 

evaluated at the equilibrium density. (Generally cW(r) is a functional of p(r).) Making 
the approximation that these identities, Eqs. ()17|) and (|18p. valid for the equilibrium fluid, 
hold also for the non-equilibrium fluid and substituting into Eq. ([15)1 . we obtain the DDFT 
equation, stated without justification by Evans (see Eqs. (166) and (167) of Ref.-), and 
derived much more convincingly by MT^, i.e. 



! 0p(r,t) = 

dt 



(19) 



S P (T,t) 

where F[p(r,t)} is the Helmholtz free energy functional: 

F[p(r,t)} = k B T J drp(r,t)[ln(p(r,t)A 3 )-l] 

+ F ex \p{v,t)\ + J drV ext (r,t)p(r,t). (20) 

The first term is the ideal gas free energy; A is the de Broglie wavelength. In obtaining Eq. 
()19|) by using Eqs. (fT7|) and which are strictly equilibrium results, we are effectively 
assuming that the correlations between the particles when the fluid is out of equilibrium are 
equivalent to those for an equilibrium fluid with the same one-body density profile p(r, t) . 
We can obtain further insight into the status of Eq. (fT^j) by re-writing it as follows: 

r- 1 ^!^ = v.[ P (r,t)vMM)], (21) 

where — V/i(r,t) = — V(5F[p]/5p) is the net driving force acting on a particle located at 
(r,t). The chemical potential obtained from Eq. (J2*U|) has three contributions: 

5F[p(r, t)] . . , N , N , . , nnS 

■j^r ^ = A*( r )*) = Vid{*,t) + Pint{*,t) + Pext{r,t). (22) 



The first contribution is the ideal gas entropic term /^(r, t) = A^Tln A 3 p(r, t), the second 
contribution, pj n t(r, t) = -^Tc^^r), is that due to the interactions with the other particles 
in the fluid, and the final term is simply the external potential, /j, ext (r,t) = V ext (r,t). As 
noted by MT 3 , for non-interacting particles, fii nt = 0, and Eq. (J21j) reduces to the exact 
equation for the diffusion of an ideal Brownian gas, i.e. 

r- 1 ^^ = k B TV 2 p(r,t) + V.[p(r,t)VVUr,t)]. (23) 

For an inhomogeneous interacting fluid in equilibrium the density profile satisfies 5 

r , Y = u = constant, (24) 

w) 

where the Helmholtz free energy functional F[p] is given by Eq. ()2()j) with p(r,t) replaced 
by p(r) and V ext (r,t) replaced by V ext (r), so there is no net driving force on the particles. 
It follows that for a time independent external potential: V ext (r,t) — > V ext (r) as t — > oo, 
regardless of the initial profile p(r, t = 0), Eq. (jl9j) will yield the same one-body density 
profile in the limit t — > oo as does the equilibrium DFT (i.e. the solution to Eq. (|24jl ) for 
the same external potential. The present derivation provides, we believe, additional insight 
to that of MT 3 ' 4 into the physics incorporated into Eq. (flH|l (the key DDFT equation). 
Assuming one has an accurate expression for the Helmholtz free energy functional Eq. (|20j). 
and in particular for the excess Helmholtz free energy functional, Eq. (flUj) should provide 
an accurate description of the dynamics of p(r, t) for a system of Brownian particles. 



IV. SPINODAL DECOMPOSITION FROM THE DDFT EQUATION 

In this section we apply the DDFT derived in the previous section to fluid spinodal de- 
composition. Since the basis for the DDFT is the Smoluchowski equation, an equation of 
motion for Brownian particles, we expect our theory to be particularly relevant to spinodal 
decomposition in colloidal fluids, rather than molecular fluids. In colloidal fluids friction 
with the solvent results in a much faster equilibration of the momentum degrees of freedom 
compared with those of the particle positions. However, since we do not explicitly include 
the particles of the solvent in which the colloids are suspended, our theory neglects the hy- 
drodynamic interactions that may be significant to spinodal decomposition in some colloidal 
fluids. 



S 



When a (colloidal) fluid which exhibits liquid-gas phase separation (or more generally 
fluid-fluid phase separation) is quenched to a state point in the region of the phase diagram 
where there is coexistence, the fluid can phase-separate in two distinct ways. The first 
mechanism is that which occurs when the state point to which the fluid is quenched is near 
to the binodal. In this case phase separation generally occurs via nucleation of droplets of 
one phase forming in the other phas o 32 i 33 . For example, if the fluid is quenched to a state 
point inside the binodal on the liquid side, then bubbles of the gas phase can nucleate in 
the bulk of the metastable liquid. 

However, if the fluid is quenched to a state point well inside the binodal, then a dif- 
ferent mechanism for phase separation is possible: spinodal decomposition. Spinodal de- 
composition is characterised by the exponential growth of density fluctuations of certain 
wavelength o 32 ! 33 . In (mean-field) theoretical descriptions, liquid-gas phase separation is de- 
termined by the occurrence of a van der Waals loop in the Helmholtz free energy per particle, 
f(v), where v is the volume per particle. Spinodal decomposition is predicted to occur in 
regions of the phase diagram where (d 2 f / '8v 2 )t < 0, i.e. regions where the isothermal com- 
pressibility xt is predicted to be negative. The boundary to this region, the spinodal, is 
defined by the locus of (d 2 f/dv 2 )T = in the phase diagram. Experimentally there is not 
necessarily a sharp distinction between regions where phase separation occurs via nucleation 
and via spinodal decomposition. However, in a deep quench far from the spinodal, spinodal 
decomposition is the mechanism generally expected for phase separation. 

In a fluid undergoing spinodal decomposition three different regimes can be distinguished. 
For early times after the quench, the amplitude of the density fluctuations are small and 
theories linear in the density fluctuations such as the well known Cahn-Hilliard theor y 34 ' 35 
(see also Refs . 29 ' 32 ! 33 ' 36 ! 37 ) provide a good description of this (early) stage of spinodal decom- 
position. At intermediate times the density fluctuations can be large, but sharp interfaces 
between domains of gas-like and liquid-like regions have still not formed 2 ^,. At later stages 
there are sharp interfaces between domains of liquid and gas and successful theoretical de- 
scriptions of the later stage dynamics of spinodal decomposition, such as the Allen-Calm 
theory 3 - 9 - (see also Refs . 32 ' 33 ' 38 ), focus on the dynamics of the interfaces. 

First, in Sec. IIV A[ we shall use the DDFT formalism to derive a generalisation of the 
(linear) Cahn-Hilliard theor y 34 ' 35 , similar to that described in Refs . 29 i 36 i 37 , for the early 
stages of spinodal decomposition, when the density fluctuations are small. Then in Sec. 
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II V Bl we will proceed to derive a non-linear theory which we believe may be applicable for 
the dynamics of spinodal decomposition at both short and intermediate time scales. Results 
of explicit calculations for a model fluid are given in Sec. II V CI 



A. Early stages of spinodal decomposition 

We consider spinodal decomposition in the bulk of a fluid, so we set V ex t(r,t) = in 
Eqs. (dnj) and (|20~j) and consider small density fluctuations p(r, t) = p(r, t) — p b about the 
bulk fluid density, pb, i.e. we are considering a homogeneous fluid which has been rapidly 
quenched to the region of the phase diagram inside the spinodal (e.g. from A to B in Fig. 

and we seek those wavenumbers k for which density fluctuations grow. From Eqs. (|18|) . 
(fUl and (1201) we obtain: 

{Tk B T)-^ d -^fi = V 2 p(r,t) - p fc V 2 c«(r,t) 

-V.[p(r,t)Vc«(r,t)]. (25) 

This approach is basically that of Refs.^i. In Ref.— Evans writes down Eq. (^HJ), and then 
linearises Eq. (|23J) in p by Taylor expanding about the bulk fluid value, giving 

,5cW(r) 



c «(r) = c«(oo) + J dr' 



5p{r') 



p(r',t) + 0(p 2 ), (26) 



Pb 

where cW(oo) = c^fpj,] = —fi^ex and p ex is the excess chemical potential. The second term 
simplifies by recalling^ 

*c«(r) _ 6 2 F ex [p] 



5p(r') 5p(r')5p(r) 



= c^(|r-r'|;p 6 ), (27) 

for a homogeneous fluid of spherically symmetric particles. For an equilibrium system 
c( 2 )(r;p&) is the Ornstein-Zernike pair direct correlation function of the fluid of density 
Pb- Substituting Eq. (|2Hj) into Eq. (|23|) . keeping only terms that are linear in the fluctuation 
p, we obtain^^: 

(r^T^^M = v 2 p(r,t) 

-p b V 2 [ J dr' C ( 2 )(|r-r'|;p fe )p(r',t)]. (28) 
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Fourier transforming yields an equation for the time evolution of the different Fourier com- 
ponents 

p(k,t) = J drexp(ik.r)p(r,i), (29) 

and one obtains 

(T kB T)^^M- = -k 2 p(k,t) + p b k 2 c(k)p(k,t), (30) 

where c(k) = J dr exp(zk.r)c^ 2 ^(r; /?&). The solution of Eq. (j%Dj) is 

p(k,t) = p(k,0)exp[R{k)t}, (31) 

where R(k) = —YksT k 2 (l — pbd(k)). For an equilibrium fluid, at a state point outside the 
spinodal, S(k) = (1 — pbc{k))~ l is the static structure factoid and, since for an equilibrium 
fluid S(k) > for all values of k, it follows that outside the spinodal R(k) < for all values 
of k. From Eq. (|31|). all Fourier components will decay implying, of course, the fluid is 
stable. 

On approaching the spinodal from the single phase region one finds that S(k = 0) — > oo; 
at the spinodal (1 — pi,c(k = 0)) = 0. Within the mean-field (van der Waals-like) theory 
of fluids to be described below, we find that inside the spinodal the quantity R{k) can have 
positive values for k < k c , where the value of k c depends upon how far into the spinodal 
region one has quenched, see for example Fig. El which employs a particular approximation, 
namely Eq. (jUJ), for c(k). k c is obtained as the solution to the equation pbc(k c ) = 1. Thus 
we find that inside the spinodal region there will be some density fluctuations with a wave 
number k < k c whose amplitude will grow exponentially^ 12 ^ 1 ^ 1 ? 41 ?^^ 1 ^ . The deeper the 
quench into the spinodal region, the larger the value of k c can be. 

This general picture of the exponential growth of certain Fourier components (modes) 
was obtained originally by Cahn and Hilliard 34,35 who derived an explicit approximation for 
the function R(k). Cahn-Hilliard theory for spinodal decomposition is usually derived by 
considering the continuity equation ((2)), together with the following approximation for the 
current 

jMH-mv^MI, < 32 > 

where M is a mobility constant and the functional F is chosen to be of the 4 Ginzburg- 
Landau form. Generalising slightly we assume that the free energy functional has the square- 
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gradient forroi: 



F sg [p(r,t)\ = / dr 



fo(p(r,t)) + -K\Vp(r,t)f 



(33) 



where fo(p) = pf{p) is the Helmholtz free energy density for the homogeneous fluid of 
density p and we treat K as a positive constant. From Eq. (J2J) one then finds 



dp(r,t) 
dt 



MV 2 



df (p(r,t)) 
dp{r, t) 



- KV 2 p(r,t) 



This equation is then linearised about the bulk density p b to obtai: 



,32,36. 



dp(r,t) 
dt 



MV 2 



dp 2 



p(r,t). 



(34) 



(35) 



On Fourier transforming Eq. (|3*5Jl one obtains dp(k,t)/dt = R sg (k)p(k,t), where R sg {k) = 
-Mk 2 (Kk 2 + (d 2 f /dp 2 ) Pb ), the solution to which is Eq. flHH), with R(k) replaced by R sg (k). 
Note that (d 2 fo/ 'dp 2 ) Pb is negative inside the spinodal. It is clear from this argument that 
Cahn-Hilliard theory can therefore be viewed as a special case of the more general linear the- 
ory presented in the earlier part of this subsection^. By employing a free energy functional 
more accurate than (|33|) one should be able to incorporate short wavelength density fluctua- 
tions as well as the long wavelength, k — > 0, fluctuations accounted for by a square-gradient 
approach. 



B. Spinodal decomposition at intermediate times 



In order to go one step beyond the lowest order (linear) theories described above, we con- 
sider an approximate excess Helmholtz free energy functional, obtained by Taylor expanding 
the excess Helmholtz free energy about the uniform density. By integrating Eq. (J26|) and 
omitting terms beyond 0(p 2 ) we obtain 



F ex [p] = F ex [p b ] + p ex J drp(r,t) 
-^JdvJ dv'~p(r,t)p(v',t)c^(\r-r'\;p b ). 



(36) 



This truncated quadratic density expansion is often used in the theory of inhomogeneous 
fluids in equilibrium 1 . Using Eq. (|36|) in Eq. (|2~5j) we find: 



(r^T)- 1 ^^ = v 2 p(r,t) 



-V. 



(p b + p(r,t)) J dr'p(r',t)Vc( 2 H|r-r'|;p b ) 



(37) 
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The first two terms on the right hand side are those entering the linear theory (|28j) while 
the third term is the only non-linear one which arises for the functional (j3T?j) . Fourier 
transforming Eq. (J37j) we obtain: 

+ J^f S dk ' k W t )Wl k - ,t 'l' t )' ( 38 ) 

which should be compared with Eq. (j3*Uj) : there an additional term on the right hand side 
that is non-linear in p(k,t). This term describes the coupling between the different Fourier 
components of the density fluctuations (modes )^. An equation almost equivalent to Eq. 
(JHBj) was derived recently by considering the mobility M to be a linear function of the order 
parameter, rather than a constant, in a Cahn-Hilliard treatment - see Eq. (10) in Ref.— ; 
we shall discuss this further in Sec. IIVDI In order to proceed further we must assume a 
particular form for the Helmholtz free energy functional, from which we can obtain c(k) and 
thus solve Eq. (J3~%j) numerically. 



C. Results for a model fluid 



We consider a fluid composed of particles interacting via pair potentials of the form 

v 2 {r) = v hs (r) + v at (r), (39) 

where 

{oo r < a 
(40) 
r > o 

is the hard-sphere pair potential, and the attractive part of the pair potential has a Yukawa 
form: 

a\ 3 exp(-Ar) 

Mr) = ; (41) 

where a and A are positive constants. Provided the decay length A -1 is sufficiently large 
this model fluid exhibits stable, with respect to freezing, liquid-gas phase separation. Pair 
potentials of this form are often used as crude models for simple fluids but they could be 
used to model the effective (depletion) potential between the colloids in a colloid-polymer 
mixture solution 43 by choosing the pair potential parameters in Eqs. ()40|) and ()41|) to mimic 
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the well-known Asakura-Oosawa potentia l 43 i 44 i 45 i 46 . In calculations for our model system we 
approximate the excess Helmholtz free energy functional by 



F ex [p] = F e h :[p] + IJdvJ dv f p(v,t)p(v',t)v at (\v-r'\), 



(42) 



where F^[p] is the reference hard-sphere Helmholtz excess free energy functional and attrac- 
tive interactions are treated in a mean-field fashion 1 . If we employ the Rosenfeld Fundamen- 
tal Measure Theor y 4 7i 48 i 49 for F^[p], this non-local functional generates the Percus-Yevick 
pair direct correlations functions in a hard sphere fluid. Thus, using Eq. ()27)1 we obtain the 
following simple (Random Phase) approximation 



.(2) 



(r;p b ) = c PY (r;p b ) - (3v at (r), 



(43) 



where cpy{r) is the Percus-Yevick (PY) approximation 31 for the hard-sphere pair direct 
correlation function. With this choice the Fourier transform of c^(r) can be carried out 
analytically: 



c(k) = c PY (k) + 



pa\ 2 
A 2 + k 2 ' 



(44) 



where cpy{k) is given 



c PY {k) = - 47ra 3 



a + 2(3 + 4 7 24 7 



Q 



sin(g) 



+ 



+ 



a + (3 + 7 2/3 + 127 24 7 



cos(g) 



24 7 
q 6 



where q = ka, and the coefficients 



a 



7 



2/3 
g 4 



[l + 2r,f 



(45) 



(1-r^) 4 ' 
-6t/(1 + t//2) 2 

(l-^) 4 ' 



(46) 



2(1 -nY 

depend upon r\ = np b a 3 /6, the packing fraction. The phase diagram for our model fluid 
is displayed in Fig. ^ The liquid-gas binodal is calculated by performing the common 
tangent construction on the Helmholtz free energy per particle for the bulk fluid, f(pb) = 
fpyiPb) — Pb^/2^ obtained from Eq. (j4^j) ; fpyiPb) is the PY compressibility approximation 
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for the hard-sphere Helmholtz free energy 3 ^. The spinodal is the locus of points for which 
d 2 f/dv 2 = 0, where v = 1/pb- For simplicity we choose to set the inverse length scale 
parameter in the attractive part of the pair potential equal to the hard-sphere diameter, i.e. 
A -1 = a. The quantity a = J drv at (r) determines the energy scale. We expect the present 
model fluid to exhibit a freezing transition for large r], but we do not consider this here. 

We are now in a position to use our approximate form for c(k), obtained from Eqs. ()44|) 
and (jUiJ), together with Eq. (JHHj) . to calculate p(k, t) during the early and intermediate times 
of spinodal decomposition. We assume that at t — p(k,t) takes a small, constant, positive 
value for all values of k and choose p(k, t = 0) = 1CT 8 . The later time dynamics is insensitive 
to the specific choice of value for p(k, t = 0) because for short times, when Eq. (jHTJ) describes 
the dynamics of spinodal decomposition, Fourier components (modes) with wave number 
k > k c do not grow and are exponentially damped, whereas components with k < k c grow 
exponentially with the initial value p(k, t = 0) as a prefactor. Thus choosing a different 
value for p(k, t = 0) is equivalent to a shift of the time axis. For a quench to the state point 
ksTa^/a = 0.05 and r] = 0.2 (i.e. the quench to point B, in Fig. the early time growth 
gives a single peak in p(k,t) at ko ~ 0.6 corresponding to where the maximum of R(k) 
in Eq. (p?Tj) occurs - see the inset to Fig. El We find that R{k) obtained using Eqs. (}4"4"|) 
and ()45|) for c(k) has a similar form to that extracted from molecular dynamics simulation 
results for a Lennard- Jones fluid deep inside the spinodal region^. Moreover, the overall 
shape of R(k) and the values of k c a are in keeping with results for the Lennard- Jones fluid 
(for similar state points) obtained by Evans and Telo Da Gama^, using a theory equivalent 
to the present, and by Abraham 37 , using a perturbation theory approach. We also plot in 
Fig. |2]the function S(k) = (1 — p b c(k))~ l ; the negative portion of this function corresponds 
to wave numbers k for which density fluctuations grow exponentially in the early stages of 
spinodal decomposition. 

Some typical plots of p(k, t) at early times are displayed in the inset to Fig. El for a quench 
to point B in Fig.^ The plots are for the reduced times t* = kBTTa 2 t = 45, 50 and 55 in the 
inset and 60, 65, 67 and 69 in the main figure. The results obtained from the linear theory, 
Eq. (|nUJ), and the non-linear theory, Eq. ()38|) . are indistinguishable for the three earliest 
times (see inset). However, at later times the linear theory (dashed line) continues to give 
just a single peak in p(k, t) that grows exponentially, whereas the non-linear theory(solid 
line), which includes the effect of coupling between Fourier components with different wave 
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numbers, shows that components with wave numbers different from that predicted by the 
linear theory can also grow. We see that the effect of the coupling incorporated into the 
non-linear theory becomes increasingly significant at intermediate times, producing first a 
shoulder which grows into a bump in p(k, t) at a larger wave number, ka ~ 0.8 than that of 
the main peak which first appears in p(k, t) at early times. 

D. Discussion 

In simulation studies of spinodal decomposition the quantity that is often measured in 
order to characterise the fluid is the structure factor 



where (...) is an ensemble average over different realisations of the stochastic noise in the 
interval up to time t. For small values of k one finds thalj^ 9 -: 



where A(t) is a (small) wave number independent baseline. Thus, to a reasonable approx- 
imation S(k,t) oc [p(k,t)} 2 2& (see also Ref.^). Simulation studies of the early stages of 
spinodal decomposition in model colloidal fluids such as that in Ref.— , where they consider 
a fluid composed of particles interacting via the Lennard- Jones pair potential, display the 
growth of a single peak in S(k,t). This is, of course, a general feature of the very early 
stages of spinodal decomposition and is found in many other system a 29 ' 32 . Thus, our results 
for short times showing the growth of a single peak in p(k,t), and therefore also in S(k,t), 
are in keeping with simulation studies. 

The development of a second peak (shoulder) in S(k,t) at intermediate times after the 
quench, at a larger value of k was observed in the two-dimensional calculations of Mao 
et al. 42 , where they considered a non-linear extension of Cahn-Hilliard theory for polymer 
mixtures. The order parameter for their theory is the deviation of the concentration c(r, t) = 
C(r, t) — C (i.e. c(r, t) replaces p(r, t) in our theory). They used the approximation that the 
mobility M = M + Mic(r, t), i.e. a linear function of c(r, t), together with a square-gradient 
Helmholtz free energy functional in Eq. By linearising 5F/5C about the average 

concentration Cq in the manner leading to Eq. (jHo^) along with the continuity equation ©, 




(47) 



S(k,t) ~ A(t) + hp(k,t)] 2 , 



(48) 
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Mao et al— obtain a non-linear theory that is very similar in structure to that which we 
obtain using Eq. (0) with a constant mobility together with the truncated functional (|36j). 
The main difference between the two approaches, other than the choice of approximation for 
the Helmholtz free energy functional^, is that Mao et al. 42 have an adjustable parameter, 
the ratio Mx/Mq, that allows a tuning of the mobility for their polymer blend, whereas in 
our approach we are effectively restricted to the choice Mi/Mq = 1/pt- Their results also 
show that the second peak (shoulder) in S(k, t) results from considering a theory that is non- 
linear (second order) in the order parameter. The other significant difference from Ref.— is 
that the present theory uses a microscopic non-local functional (|3fi|) . and therefore includes 
the effects of interparticle correlations more accurately than a gradient (Ginzburg-Landau) 
theory such as that used by Mao et alM. Nevertheless, in the early and intermediate times 
of spinodal decomposition, where sharp interfaces between gas-like and liquid-like domains 
have not yet formed, we should expect good qualitative agreement between our results and 
those from the non-linear Cahn-Hilliard gradient theory of Mao et alM. Mao et al. also found 
that the presence of the non-linear terms in their theory resulted in a significant change in 
the connectivity of contour plots of the order parameter from that found in the absence of 
these terms (Mi = 0) - see Figs. 3, 5 of Ref.—. On the basis of our present analysis we would 
expect this observation to be a general feature of spinodal decomposition at intermediate 
time scales. 



V. CONCLUDING REMARKS 



This paper falls broadly into two parts. In the first part we provided an alternative 
derivation of the DDFT developed by MT 3 4 . Our derivation elucidates the physical approx- 
imations made in order to construct the theory. The starting point is the Smoluchowski 
equation (JHJ), the (generalised) Fokker-Planck equation for the probability distribution func- 
tion P(r N ,t) corresponding to the Langevin equation of motion, Eq. (fTUj) . This stochastic 
basis for the theory means that the theory should be applicable to colloidal fluids where 
because of interactions with the solvent particles the momentum degrees of freedom of the 
colloids equilibrate much faster than the positional degrees of freedom. In atomic fluids the 
equilibration timescale for the momentum degrees of freedom can be of the same order as 
that for the relaxation of positional degrees of freedom and therefore the present theory may 
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break down for atomic fluids. However, because the correct equilibrium limit is built into 
this theory, i.e. when dp(r,t)/dt = the theory is equivalent to Eq. (|24p. it is feasible that 
the present theory would give reasonable results for atomic fluids, assuming that one also 
has a reliable approximation for the excess Helmholtz free energy functional, F ex [p\. The 
fact that the correct equilibrium limit is built in is, we believe, one of the most appealing 
features of the theory. 

The further approximation (beyond assuming that the Smoluchowski equation holds) 
made in deriving the DDFT, Eq. (fT§|). is to assume Eqs. (JT7J) and (JTHjl . exact equilibrium 
results, are also valid for the non-equilibrium fluid. This approximation is equivalent to 
assuming that the two-particle, three-particle and higher order correlations in the non- 
equilibrium fluid are equivalent to those in an equilibrium fluid with the specified one-body 
density profile. We expect that this approximation is a reasonable one to make and will not 
result in significant failures of the theory, as long as other approximations can be justified. 
We therefore believe that Eq. (|19|) may well provide a reliable account of the dynamics for 
a variety of different fluids, provided of course, that one has an accurate approximation for 
the appropriate F ex [p\. 

The second part of the present paper is concerned with the application of the DDFT to 
the problem of spinodal decomposition. Our key result is Eq. (|38|). For early times after 
quenching the fluid, when density fluctuations are small in amplitude, our theory reduces 
to the linear theory of Evans and Telo Da Gama^. Furthermore, this linear theory reduces 
to the Cahn-Hilliard theory if we were to use a square gradient Helmholtz free energy 
functional Eq. (|3*3*|) rather than the more accurate non-local functional, Eq. ()36j) . Cahn- 
Hilliard theor y 32 i 33 i 34 i 35 is known to provide a successful qualitative description of the early 
stages of spinodal decomposition. Since the present theory incorporates all the effects that 
Cahn-Hilliard theory describes, and provides a more accurate treatment of short wavelength 
correlations 36 it should be a reliable quantitative theory for the early stages of spinodal 
decomposition. 

However, the key feature of our theory is that it incorporates a "mode- coupling" term 
(final term on the right hand side of Eq. ([38))) which describes the coupling between different 
density fluctuation modes, an effect which becomes important in the intermediate stages of 
spinodal decomposition. The results of calculations including our "mode- coupling" term 
are in qualitative agreement with those from a recent non-linear extension to Cahn-Hilliard 
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theory— We believe that further work needs to be done in testing the predictions of our 
theory. In particular, a comparison with Brownian dynamics simulation results for spinodal 
decomposition would be very useful; we do not know of relevant simulations which go beyond 
the very early stages of spinodal decomposition. Of course, running simulations that get 
into the stage of spinodal decomposition for which our results exhibit deviations from those 
of the linear theory could be computationally expensive. 

We conclude by mentioning one important conceptual issue regarding the input to our 
theory of spinodal decomposition. Clearly Eqs. (|5U|) or (|3H|) require c(k) as input. As em- 
phasised earlier, for an equilibrium state c(k) is simply the Fourier transform of the pair 
direct correlation function and, as such, it can be obtained from integral equation theories^ 
or from simulations. In Ref.— results for c(k) calculated using Percus-Yevick theory for a 
Lennard- Jones fluid were extrapolated into the unstable spinodal region 54 ' 55 . Such a proce- 
dure is fraught with uncertainty and is difficult to justify. The present DFT formulation 
of the theory avoids such problems, is defined as the second functional derivative of 
— /3F ex [p] and provided one makes a division of -F ex [p] into a repulsive reference contribu- 
tion plus an attractive (mean-field-like) contribution, as in Eq. (|42jl . there is no difficulty 
in calculating c^(r;p&) inside the spinodal. We do not need to make any extrapolation. 
One could envisage treating attractive interactions in a more sophisticated fashion 1 - than in 
Eq. (|4~2*j) but provided the basic division is maintained one should expect to obtain similar 
results for c(k). 
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FIG. 1: Phase diagram for a fluid composed of particles interacting via pair potentials of the 
form in Eqs. (|39B41j) . calculated from the free-energy of Eq. (|42|) . rj = irpbcr 3 /6 is the packing 
fraction and ksTa 3 /a is the (reduced) temperature. The path from A to B denotes the quench 
corresponding to the results in Figs. [21 and [U 
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FIG. 2: The function S(k) = (1 — pbd(k))^ 1 calculated using Eq. (|H|) at the state point with 
r\ = 0.2 and ksTa 3 /a = 0.05, corresponding to point B in Fig. ^ Note that for k < k c , where 
k c a ~ 0.8, S(k) < 0. In the inset we plot R(k)a 2 (3 /V = —k 2 a 2 /S(k), the factor appearing in the 
exponential in Eq. (jHlj) . In the initial stages of spinodal decomposition density fluctuations with 
wave numbers k < k c grow exponentially, whereas for k > k c the fluctuations are damped. 
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FIG. 3: Plot of p(k,t) for a quench to the state point ksTa^/a = 0.05 and r\ = 0.2, point B in 
Fig. □ p(k, t) is shown for times t* = k B TTa 2 t = 45, 50, 55 (in the inset) and t* = 60, 65, 67, 69 (in 
the main figure). The results obtained from the linear theory (Eq. (jSOl)) are denoted by a dashed 
line, and those from the non-linear theory, (Eq. Q38|)) by a solid line. The effect of the coupling 
between Fourier components (modes) with different wave numbers, described by the non-linear 
theory, becomes increasingly significant at later times, whereas for earlier times (see inset), the 
results from the two theories are indistinguishable. 
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